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1. MATHEMATICAL DEVELOPMENT OF HISS 


1 .0 Introduction 

NASA is active in research on supersonic combustion [1-3]. Of 
particular interest are methods for injecting hydrogen in a supersonic 
air stream in a manner which achieves effective mixing and a minimum 
pressure loss. Two arrangements of interest are shown in Figures 
(1.0-1) and (1.0-2). 

This report provides documentation of a finite-difference computer 
program called HISS (^drogen injection of a Supersonic Stream) which 
was developed by FMTS under contract to NASA for the exclusive use of 
NASA; it can be used to predict the flow properties for cases depicted 
in Figures (1.0-1) and (1.0-2). 



Figure I . Or- I 


Normal 


Injection 




Figure 1.0 — 2 Parallel Injection 


The remaining part of this chapter gives the mathematical equations 
which are used to model the flow process. The "art" involved in using 
these equations to develop a finite-difference algorithm for these 
equations is also described. 

1.1 Capabilities and limitations of HISS 

HISS is capable of calculating three-dimensional, boundary layer 
flows which are either external or internal. The geometry considered 
is flow in a fixed rectangular parallelepiped with "zero- flux" boundary 
conditions on each side wall, either a free, symmetry or solid wall at 
the top surface and either a symmetry or wall boundary condition at the 
bottom surface. It can handle injection of hydrogen ranging from 
normal to parallel to the main air stream. Five equilibrium reactions 
are allowed. The main air flow can be either subsonic or supersonic; 
however, the free stream pressure gradient is considered to be zero. 

The flow is considered turbulent and the viscosity is calculated 
by way of the "k-e" turbulence model described in [4J. Density is 
calculated via the ideal gas law. 

1 .2 Numerical method 

The numerical method used is described in detail in I5H. A brief 
discussion will be made to summarize the metFlod. 
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1.2-1 The differential equation 


The general differential equation describing the transport of a 
fluid property <|) is expressed in Cartesian coordinates as follows for 
the fixed volume under consideration: 


^ (pu^) + ^ (pv4.) + ^ (pwd>) = S4. + ^ (r^ 1^) + ^ (r* |i) 

( 1 . 2 - 1 ) 

where x, y, and z are Cartesian coordinates as illustrated in Figure (1.2-1). 
(|) is a general dependent variable. 

The governing equations which must be solved are: continuity, 
momentum equations for each coordinate direction, the energy equation, 
an equation expressing conservation of total hydrogen and two differential 
equations for turbulence properties. These equations allow calculations 
of eight dependent variables, namely u, v, w, p, f, R, k and e. All of 
these equations except continuity have the forro of equation (1.2-1) and 
the respective source terms and transport coefficients are given in 
Table (1.2-1). The calculation of pressure is described in Section 1.5. 



Figure 1.2-1 Coordinate Systems 
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For these variables, the appropriate transport and source terms 
are given in Table (1.2-1). 
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Table (1.2-1) 





















Values for the laminar Prandtl number, a, and the turbulent Prandtl 
number ^ along with C-j and Cg required in the source term of the k and 
e equations are given in Table (1.2-2). Also given are 


<l> 

"’eff,.|. 

a 

k 

1.0 

.7 

e 

1.3 

.7 

h 

.9 

.7 

f 

.9 

.7 


K 

E 

Cd 

Cl 

C2 

.42 

9.0 

.09 

1.44 

1.92 


Table (1.2-2) 

K and E which are constants required in the law-of-the-wall formulation 
described in section 1.6-4 and the constant Cj^ used to obtain the 
dissipation length scale in equation 1.3-2. 

1.3 Auxi 1 i ary i nf ormati on 

In addition to the partial differential equations the complete 
specification of the problem requires provision of auxiliary information 
of three types: boundary conditions, physical hypotheses which permit 

the calculation of diffusion coefficients as well as sources and sinks 
of each variable, and certain relationships among the thermodynamic 
and transport properties required to describe the flow. 

1.3-1 Boundary conditions 

For each of the seven variables listed in Table (1.2-li, values 
must be specified in the z = 0 plane. In addition, boundary conditions 
along th.e N and S surface Isee Figure 0.2-lD'nitfst be specified for each 
of these variables. These boundary conditions can be specified 
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as the value of the variable <(. or the flux of ^ through the surface. A 
detailed discussion of boundary conditions is given in Section 1.6. 

1.3-2 Physical hypotheses 

The gas density is calculated by the ideal gas law. Details are 
given in Section 1.7 since special procedures are required to treat mixed 
subsonic-supersonic flow. 

The laminar viscosity is calculated by the formulas presented on 
page 58 of reference £6J and will not be repeated here. 

The turbulent viscosity is calculated by solving conservation equa- 
tions for the kinetic energy of turbulent fluctuating motion, k, and the 
energy dissipation rate, e. The turbulent viscosity is then given by: 

Pk^)^ (1.3-1) 

where is a dissipation length scale which can be determined from the 
equation for e, i .e. : 

C„ k3/2 

)l = (1.3-2) 


Equations (1.3-1) and (1.3-2) are the basis by which the effective turbu- 
lent viscosity is calculated by the k-e model. 

1.3-3 Thermodynamic and transport relationships 

In order to calculate properties such as laminar viscosity or 
temperature, it is necessary to know the mass fraction of each chemical 
specie present at a given location. Thermodynamic equilibrium is assumed 
so that from the local temperature, pressji^re, and element fractions the 
specie mass fraction distribution is calculated. This technique was 
developed in this work and is reported in detail in Appendix A. 
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Another important relationship is between the stagnation enthalpy 
and temperature. This relationship is: 


R = 


En.f,h« 


U'^ + 


+ c. 


- ■'ref> 


where Cp is the mixture specific heat and is the reference tempera- 
ture for which the specie enthalpy of formation, h°^, is based. For this 
work T^g^ = 0°K. Note that Cp is defined by: 


T-T, 


ref 


/ 


h-h 

S'"' = T:T 


ref 


'ref 


ref 


where his the sensible enthalpy at the temperature T. Therefore, the 
definition takes into account the variation of specific heat with tempera- 
ture. 

The density is assumed to vary according to the ideal gas law, 
namely: 


p 


£W 

RT 


1.4 The discretization procedure 

The finite-difference equivalents of the differential equations are 
obtained by integrating the latter over the control volumes which surround 
the nodes of a grid system. For purposes of this integration, the de- 
pendent variables are presumed to vary in a prescribed manner between grid 
nodes; details of this are available in Patankar and Spalding [ 5 ]. 

1.4-1 The grid system 

The grid system used is indicated in Figure (1.4-1). It consists of 
a system of orthogonal intersecting grid lines in the x-y plane corre- 
sponding to a constant value of z. 
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Figure 1.4-1 Grid System 


J.4-2 Storage locations 


Figure (1.4-1) also indicates the fact that the u and v velocity 
components are located at points in the grid regularly displaced from 
the points where all other variables are stored. The boomerang-shaped 
envelopes enclose the triads of points denoted by the single letter N, S, 
E, W or P, and represent a unique computer storage location. 


1.4-3 Control volumes 

The control volume surrounding each grid node P, is indicated by 
dashed lines in Figure (1.4-2), and is termed the main control volume. 
Control volumes appropriate to the u and v velocity components are "stag- 
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Figure 1.4-2 Main Control Volume 


gered" from such main control volumes; the faces of the former being 
arranged to pass through grid nodes where pressures are stored. 

The control volumes corresponding to the near-boundary velocities, 
V in the case of N and S boundaries and u in the case of E and W bounda- 
ries, are arranged to be somewhat larger than their size in the rest of 
the calculation domain. Figure (1.4-3) illustrates this point. The 
near-boundary main control volumes remain unchanged. 
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U- Control Volumes 



Figure 1.4 — 3 Near Boundary Velocity Control Volume 


1.4-4 The general discretized equation 

The general discretized equation is obtained by integrating 
equation (1.2-1) over the control volume surrounding each grid node. 
Such an integration is performed after making assumptions about the 
manner in which the variable <() is distributed between grid nodes. All 




I If 

# 


dependent variables are presumed so to vary, linearly in the x and 
y directions and in castellated fashion along the z direction. The 
result of the above-mentioned operation is an algebraic equation for 
each grid location. This algebraic equation now represents the dis- 
cretized form of the balance equation for over the control volume 
surrounding the grid location. The balance equation is expressed 
simply as: 


[^D - ''u ‘^P,u] ^p) ■ ^s ^^p h] 

Ij-g ^ *^P^ " ^"^P ~ ^P 

* (*n - ♦?> - <*P - ♦s’] * ^*E - *P> - ^ <*P - ‘u’] 

( 1 . 4 - 1 ) 

where the L's represent convective contributions, T's the diffusive 
and S the source (and/or sink) contribution to the balance of (|>. The 
location of points e, w, n, s, E, W, N, S is given in Figure 1.4-2. 
These coefficients are defined as follows: 


fu ■ 

Hp,u 

Ay 

(a) 

■ 

t] n "T 


(b). 

■ 

Hs ^ 

- 

(c) 


fu - 2 '■I* 

2L^- 2L^ 21^ 

(d) 

n - 

Ay AZ 

-p 

jX _ Ay AZ 

(e) 

()i,e fix 
^ e 

w ,w fix^ 

= 

AX AZ 

p 

T-y „ AX AZ 

(f) 

S 

<l> , S fiyg 

T = r 

n <fi,n fiy^ 
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H AX Ay AZ 


(g) 


AXg = Xp - X^ 

Xy - Xp 

(h) 


Ays = Yp - ys 

Ayn = Yn - yp 

(i) 

(1.4-2) 

AX = |Xg - X^l 


(j) 


Ay = ly„ - ysl 


(k) 



where the subscripts e, w, n, s, E, W, N, S indicate location defined 
in Figure 1.4-2, Az is the forward step size and S'** is always expressed 
in linearized form as: 

Sj = + Sp<^p (1.4-3) 

On re-arranging the terms of equation (1.4-1), we have: 

<t>p = + A3'}>5 + + Ay((>y + B (1.4-4) 

where the coefficients A^l are defined as follows: 


\ - (T, - Lej/'V 

(a) 


(b) 
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(c) 


■ (^s * 

(d) 

(1.4-5) 

B = [(S|j + Sp^p) + Fy ^p^^j /Ap 

(e) 


a; ■ Tj F F . T* - L» + * l;; + F„ - Sp 

(f) 



It is possible for the ratio of convective contribution to 

the coefficient to become large on occasions resulting in the 

coefficient becoming negative. Were equation (1.4-4) to be solved 

with negative coefficients, physically implausible results would ensue. 

To overcome such a possibility, a hybrid scheme, well-known as the 

high- lateral -flux-modification (see Reference j]5j) is introduced. This 

scheme consists of replacing all the coefficients of the form T^jJ by 

T*^ as follows: 
m 


T" = 
m 


m 




(1.4-6) 


where, lL| signifies the modulus of L. 


1.5 The computational algorithm 

The computational algorithm embodied in HISS, is the SIMPLE 
(for Semi-Implicit Method for Pressure-Linked Equations) scheme reported 
in reference [5]. In brief, this scheme is described as follows. 

Approximate forms of discretized momentum equations are first 
solved with a guessed-at pressure-field. The resulting "starred" velocity 
field is used in conjunction with the discretized continuity equation to 
arrive at a distribution of "pressure-corrections ," p ' . These pressure- 
corrections are then used in correcting the pressure and velocity fields 
simultaneously. The corrected velocity fields are subsequently used in 
solving the discretized forms of all other equations. 

* Reference C53 describes an incompressible solution. The 

algorithm used in HISS has been modified to handle compressible 13 

flow as described in Section 1.7. 


Solutions to the algebraic equations themselves are obtained by the 
use of standard tri -diagonal matrix algorithm (TDMA). 

1.6 Boundary conditions 

1.6- 1 General policy 

The HISS program requires the specification of Boundary conditions 
on the N and S boundaries. The boundaries on each side are fixed as 
"zero flux" boundaries. A clear distinction is made, for all dependent 
variables, between boundary values and values internal to the domain. 

The main machinery of the program leaves the boundary values, unchanged, 
although it uses them in determining the internal values. Thus, the 
procedure is so structured that it nominally solves the fixed-boundary- 
value problem. When boundary values are not known, however, appropriate 
modifications are devised which permit the single structure to be used. 

The following sections describe such modifications. 

In general, boundary condition information can be supplied to the 
numerical calculation procedure in one of four ways. The boundary values 
of the dependent variables themselves can be modified; or the values of r 
at the boundary nodes can be suitably adjusted. Alternatively, the source 
terms for the near-boundary control volumes or the finite-difference 
coefficients themselves can be suitably modified. HISS is equipped with 
source- term modification practices. 

1.6- 2 Boundary conditions for transport equations 

Fixed Boundary values 

This is the simplest case and the user is merely required to supply 
the correct boundary values at the start of the calculation procedure; or 
else, if the 'fixed' boundary conditions vary in any arbitrarily specified 
manner along the z direction, account can be taken of such variations by 
suitably updating the boundary value of the appropriate variable. 

Fixed gradients at boundaries 

When the boundary value of a dependent variable is unknown but the 
normal gradient or the flux of the variable at the boundary is given, 
the following practice is adopted. The coefficient of the finite-difference 
equations connecting the boundary node to an internal node is arranged to 



be zero. The boundary flux is then supplied through the source term 
for the internal node. A common example of this is in the treatment of 
planes of symmetry. The only action taken here is to set the appro- 
priate finite-difference coefficient to zero. 

Retrieval of unknown boundary values 

In situations where boundary values of dependent variables are 
unknown and the boundary coefficient has been set to zero, the value 
of the dependent variable in the computer register corresponding to the 
boundary node has no significance- (The wall condition is therv represented 
by the value at the node nearest the wall!. In th.is program no action is 
taken to modify tt so that the initial value of the dependent variable will 
continue to prevail in such registers. 

1.6- 3 Boundary conditions for the pressure-correction equation 

Unlike the finite- difference equations for a general variable (p, 
the pressure correction equation is programmed with the presumption that 
all gradients of p' normal to boundaries are zero at the boundaries 
themselves. The pressure-difference coefficients D*^, and d'^ are 
initialized to zero over the entire calculation domain and remain at 
this value for the boundary nodes. 

1.6- 4 Wall functions 

The expressions for r appropriate to turbulent flow are not strictly 
valid in the vicinity of wall boundaries to the flow, where laminar vis- 
cosity plays an important role. If the near-boundary grid nodes are suf- 
ficiently far away from walls, the turbulent viscosity can continue to be 
used for internal grid nodes. However, means must be provided for the 
calculation of the correct shear stresses as well as fluxes of other 
dependent variables at the wall boundaries. Provisions for such calcu- 
lations are made in the program, and use the so-called wall-function 
concept. In this concept, the flux of a variable <|> at a wall boundary, 
is expressed as: 


’’wall 


= r . 


,wal 1 


*near wall " '^wall 
■^wall 


( 1 . 6 - 1 ) 
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where denotes the normal distance from the wall to the near-wall 

point. Values of are obtained from the presumption that in 

the region adjacent to wall boundaries, the dependent variables obey, 
for turbulent flow, a modified form of the semi-logarithmic law-of- 
the-wall. The formulae used to calculate ^^11 are provided in 
Table (1.6-1). 


Velocity components 
normal to the wall 


Velocity components 
parallel to the wall 


y"^ > 11.5 : 


+ 

y 


iln 

K 

EEy'’] 


11.5 : 

y 



0 


y"^ > 11.5 : ^ 
a. 


All other i>'s 


eff,i.(iln [Ey^] + P } 


< 11.5 : 


Table (1.6-1) 







The definition of in the table, is a generalization due to 
Spalding [5] of the conventional form, in that : 


+ 

y 





6 


P 


( 1 . 6 - 2 ) 


where 6 denotes the distance from the wall, at a location with which 
the values of p and k are associated. The constants k and E in Table 
1.2-2 are obtained from the conventional form of the law-of-the-wall : 

ti' = ^ In EEy"^] (1.6-3) 

and are given values of 0.42 and 9.0 respectively. 

The boundary conditions for k and e are provided as follows: 
the diffusion of kinetic energy k, to the wall is known to be negligible 
and is set to zero and a balance equation for k, regular in other respects 
is solved for control volumes adjacent to wall boundaries. The diffusion 
of dissipation rate s to such boundaries is more difficult to express. 
Instead of attempting to calculate , use is made of the knowledge 

that the length scale 1 varies linearly with distance from the wall, in 
the neighborhood of the wall. The dissipation rate is then calculated 
from this length scale from: 


^near wall 


k3/2 


r 3/4 
D K6 


Cl. 6-4) 


The practice adopted is to fix to the above value, without 

disturbing the general calculation procedure, in the following manner: 

^near wall 

(1.6-5) 

Sp = -L (,b) 


where L is a large number, say lO^®. 
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n 


The expression for for a general dependent variable <t>, 
is based on the expression indicated in Table (1.6-1), the value of 
P. being calculated from: 

9 


9.24 


\ 2 



*^eff,<|. 


‘^eff,<|) 


( 1 . 6 - 6 ) 


where a, and c!^. . denote the laminar and tubxilent values of 

(J) t jCp 

Prandtl/Schmidt number appropriate to the transport of (j>. 
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1.7 


The treatment of compressibility 

The special compressible-flow features of the computational 
procedure in HISS are as follows: 

1.7-1 The calculation and use of density 

Two distinct densities are computed and used at each step. The 
downstream density (pp) is used only in computing the mass flux in the z- 
direction at the downstream face of the cell; the upstream density 
(pp ij) is used for calculating lateral mass fluxes, and the upstream 
axial mass flux. The two densities are calculated as follows: 


_ 


P, U 


R T, 


P,u 


_ Pp.U '^P.UU 


Pp,U R T, 


P,UU 


(1.7-1) 


(1.7-2) 


where Wp is the local mixture molecular weight, R is the universal gas 
constant, and subscript ' UU' refers to the plane two steps upstream 
of the plane in question. 

1.7-2 The momentum equation in the z-direction 

The source term employed in the w equation depends 
on whether the local flow is subsonic or supersonic. When the flow 
is subsonic it is necessary , in order to render the equations para- 
bolic, to write the source term as a mean pressure gradient, i.e.. 


s = dF- P - Pu 
w d z sz 


(1.7-3) 


where p denotes the mean pressure which is determined for a confined 
flow from the requirements of overall continuity, as described in 
reference [5]; for an unconfined flow, F is simply the free-stream pressure. 


* The need for this practice is discussed in reference [5]. 
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For supersonic flow, in contrast, it is possible, without destroying 
the parabolic nature of the equations, to employ the local pressure to 
calculate the pressure gradient, i.e.. 


s = M = 

W 9z 6z 


(1.7-4) 


This practice, it should be observed, allows full account to be taken of 
pressure waves in supersonic flow. 

1.7-3 The pressure-correction equation 

The incompressible- flow form of the pressure-correction equation 
given in reference [5]* must, in order to handle compressible flow, be 
modified as described below. 

Account must now be taken of the effect of a pressure change at P 
on the mass flux at the downstream face of the cell. For supersonic flow 
the mass-flux change (pw)p is related to the pressure correction Pp by: 


(pw)^ = [p^ °p ■"(^)p '*'p*^ * Pp 


(1.7-5) 


where, pp* and Wp* are computed from the guessed pressure Pp*, (dp/dp)p is, 
from equation (1.7-1): 


[dp p - RTp^^ ’ 


(1.7-6) 


and 



(1.7-7) 


and is deduced from the finite-difference form of the momentum equations, 
in the manner described in reference [5]. 


* Note that in reference [5] the forward-marching direction is the x 
direction, whereas in HISS it is the z direction. 
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For subsonic flow, equation (1.7-5) must be modified to account for 
the fact that Wp no longer depends on Ppj the modified expression is: 

(Pw)p = (^]p • V • Pp‘ (1-7-8) 

The resulting equation for p' is of the form given in reference [5] 
with the coefficients modified to account for the above-described 
influences of Pn' on (pw) 

K p. 
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2. GENERAL FEATURES OF HISS 


2.0 Introduction 

In this chapter an overview of HISS is given with a description 
of the program flow and control, and definitions of important variables. 

2. 1 General structure 

The general structure of the HISS computer program is shown in 
Figure (2.1-1). The divided, vertical box on the left is the main pro- 
gram, HISS, which has the chief function of controlling the calling 
sequence of the subprograms. The starting point for HISS is at the top 
of this box. Points where call of subprograms are made in HISS are 
denoted by horizontal lines. Arrows indicate the ‘flow" of the computer 
program to and from the subprograms which are shown in boxes. The major 
calculations and logical decisions made by HISS are indicated in the 
spaces between the horizontal lines marking calling points in HISS. 

The subprograms are contained in subroutines which are not indi- 
vidually called. These subroutines can be classified in three categories: 
problem dependent, physical property dependent and invariant. BLOCK DATA 
and ALLMOD are the problem dependent subroutines while AUX is for physical 
property specification. SOLVE, STRIDA, STRIDB and PRINT are the invariant 
portions of the program. 

The remaining part of this chapter gives a general overview of the 
computer program while the next chapter details the operation of each 
subroutine. 

2.2 Fortran equivalents of main variables 

A few of the FORTRAN names used in HISS which will be required in 
the following sections, are introduced here. The dependent variables <() 
are stored in an array F, which it is convenient to consider as a three- 
dimensional array F (I,J,NV). Here I and J denote the location (respec- 
tively along the x and y directions) and NV identifies a particular 
variable. The three velocity components and the pressure-correction are 
included in the F array; however, for ease in understanding, separate 
arrays U,V,W and PP are also used and made equivalent to parts of F as 
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START 

HISS 



Figure 2.1-1. General Structure of HISS Computer Program 
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follows: 


U(I.J) 

= FCI.J.NVU) 

V(I.J) 

£ Fd.J.NVV) 

W(I,J) 

- F(I,J,NVW) 

PP(I,J) 

= F(I,J,NPP) 


The identifiers NVU, NVV, NVW and NPP, also used to identify U,V,W and 
PP elsewhere in the program, are assigned values 2,3,4 and 1. The largest 
values of I,J and NV for which storage is provided in the program are 
denoted by IMAX, JMAX and NNV respectively and assigned values 12, 34 and 

9. 

Although two- and three-dimensional arrays have been mentioned above, 
the computer program formally uses one-dimensional arrays, whose subscripts 
are calculated, each time they are required. Thus, F(I,J,NV) is referred 
to as F(IJNV), where: 


IJNV = I + JM(J) + NFM(NV) 


( 2 . 2 - 2 ) 


the arrays JM and NFM being calculated once-and-for-all from: 


JM(J) = (J - 1) * IMAX 
NFM(NV) = (NV - 1) * IF4AX * JMAX 


(2.2-3) 


Also W(I,J) is referred to as W(IJ) where: 


IJ = I + JM(J) 


(2.2-4) 


The four neighboring points of the location IJ, corresponding to the 
points of the compass are referred to as IJj^, IJ^, IJE^ and IJ^. These 
are calculated as: 
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(2.2-5) 


UN = 

U + 

I MAX 

US = 

U - 

I MAX 

UE = 

U + 

1 

UW = 

U - 

1 


It is easy to see that points further removed become: 
= UN + 1 

U^ = US - 1 etc. 


For a given problem, all members of the F array for which provision is 
available, may not be required to be solved for. Furthermore, some of 
these variables may be obtained from algebraic equations and not from 
the solutions to the partial differential equation. To provide for 
these alternatives, use is made of an array ISOLVE(NV). For ISOLVE(NV) 
equal to zero, the differential equation is not solved; solution is 
obtained for values of ISOLVE(NV) greater than zero. It is left to the 
user, to make further use of this facility. 

Other arrays directly related to members of the F array are: 

IPRINT(NV), FU(U,NV), TITLE(K,NV), FLUXN(I,NV), FLUXS(I,NV). Values 
of F(IJ ,NV) are printed out if IPRINT(NV) is equal to 1; otherwise, a 
printout is not obtained, FU(U,NV) stores the upstream value of 
F (U ,NV) during the iteration cycle at any given integration station. 
FLUXN(I,NV), stores the fluxes of the variable NV, on the I^orth boundary 
and FLUXS(I,NV) is for the ^outh boundary. TITLE (K,NV) with the integer 
K taking values 1 to 9 , stores a 36 character alphanumeric title of the 
corresponding variable NV. This is used to identify the values of the 
variables in the print-out. 

The quantities k, e, h, f, and T form the remaining variables of the 
F family. They are identified by the indices NVK, NVD, NVH, NVf, and NVT 
respectively and occupy NV locations 5, 6, 7, 8 and 9 in the array F(U,NV). 

No equation is solved for T; however, it is convenient to store it at F(U, 9). 
Quantities p, p and r are stored as P(U), RHO(U) and GAM(U). They are 
not members of the F array, but are included along with F in a COMMON 
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statement, and immediately following it. They can therefore be referred 
to as the (NNV + 1)^^, (NNV + 2)^^ and (NNV + 3)^^ members of the F array, 
as indeed they are in subroutine PRINT. The integers NVP, NRO and NMD 
are used to identify them. 

Although, in general, the values of r are different for each dependent 

variable, storage provision for only one set of r's is retained. 

Hence, at a given stage in the program, the GAM array contains the 
values of r, appropriate to the dependent variable under consideration. 
Similarly, provision is made for the retention of one set of coeffi- 
cients of the finite-difference equation at a given stage. Thus, AXP(IJ), 
AXM(IJ), AYP(IJ), AYM(IJ) and AZ(IJ) respectively represent the coeffi- 
cients A^, A^, A|Jj, A^ and Ly in equation (1.4-5). Similarly, SU(IJ) 

and SP(IJ) represent S,, and S„ respectively in equation (1.4-3). 

U^. ^ 

The quantities like Dp in reference £5] are stored as DU(IJ), there 
being similar arrays, DV(IJ) and DW(IJ), for the corresponding quantities 
associated with the v- and w- momentum equations. 

KBCN and KBCS denote the type of North and South boundary respec- 
tively as follows: 1 indicates symmetry, 2 indicates solid wall, and 
3 indicates a free boundary. 

2.3 Program control 

Program control, including start, internal monitoring and stop 
functions, is achieved through program DAVE. 

Using the information supplied through BLOCK DATA, the first part 
of DAVE calculates and assembles all the information about the grid system 
through calls to STRID0 and STRIDl. Following this, in the same part, a 
call to BEGIN, supplies the initial conditions for all the dependent 
vari ables. 

The second part begins with the statement: 60 CONTINUE. Following 
this, mass sources, step length and the location of the station for the 
new calculation, are determined. This is followed by calls to SPECIE, 
DENSTY and UPSTRM. SPECIE calculates the chemical specie composition 
from element compositions, temperature and pressure. DENSTY is a member 
of the physical modelling subroutine AUX, wherein the fluid density p, 
is calculated. UPSTRM stores the upstream variables required by the 
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computation in the array FU. After calculating the overall balance of 
the major F arrays a call is made to VISCOS which calculates r. If 
ISTEP equals NPJUMP profiles of variables specified by IPRINT(NV) = 1 
will be printed out. Otherwise, only overall balance information and 
station location will be printed out. 

The next section is for the calculation of the velocity and pressure 
fields. First, a check is made to see if injection occurs. If this 
test is positive, a call is made to INJMOD which provides the appro- 
priate boundary conditions through source terms. Next a check is made 
to determine whether the station for calculation is immediately down- 
stream of injection. If this test is positive the number of sweeps on 
the pressure correction equation and the step length are modified to 
increase accuracy and stability of the program. The calculation is 
stopped at this point if ISTEP > LASTEP. Next a check is made to deter- 
mine if ISTEP = 0 which if positive, results in iteration on the 

hydrodynamic and turbulence equation to improve the starting profile. 

After this series of tests, a call is made to STRID3 which yields 
the velocity and pressure fields. A call to STRID4 performs similar 

calculations for all other dependent variables. If ZU < ZLAST, control 

is returned to a point in DAVE just after the CALL BEGIN statement 
and the process repeated for the next station. 

2.4 Detailed list of program variables 

All variables used in common statements are defined in Appendix B. 
Variables which appear locally are defined by their use in the computer 
program and are not defined in this report. 
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3. INVARIANT PORTIONS OF HISS 


3.0 Introduction 

Many of the calculations do not need changing for different 
boundary conditions, etc. In this chapter the function of each in- 
variant portion of HISS is given. Again important variables are 
defi ned. 

3.1 Calculation of grid quantities 

The two sets of grid quantities specified by the user in BLOCK 
DATA are LCV, MCV and ZETA(I), AGEOM. The former represent the number 
of (f>-control volume faces in the x and y directions and the latter the 
non-dimensional coordinate distributions of the .^-locations in the x 
direction. 

Given the above information, STRID0, the first member subroutine 
of STRIDA, computes the maximum number of grid nodes in the I and J 
directions. This is done in the following sequence: 


L = LCV + 1 
M = MCV + 1 
LPl = L + 1 
MPl = M + 1 


(3.1-1) 


It is emphasized here that users must ensure that LPl and MPl are always 
less than or equal to IMAX and JMAX respectively. The latter represent 
the maximum dimensions of all variables in the respective directions and 
are given values accordingly in BLOCK DATA. STRID0 is concluded with the 
specification of the integer arrays JM and NFM in accordance with equation 
(2.2-3). 

The second member of STRIDA is STRIDl. STRIDl is called once from 
HISS. It is in STRIDl that the physical coordinates x and y are computed 
from the values of AGEOM and ZETA, as follows: 
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X(I) = ZETA(I) * BXE 

Y(J) = Y(J - 1) + AGEOM ** (J - 2) * DELY 


(3.1-2) 



where DELY = BYN*(1.0 - AGE0M)/(1,0 - AGEOM ** M). Here BXE and BYN 
represent the width and height of the calculation domain respectively; 

as such, they are to be specified by the user in BLOCK DATA. The above 
specification of Y(J) makes use of a geometric series for specifying 
successive intervals between Y(J) and Y(J - 1). Larger values of AGEOM 
result in the grid being "crowded" closer to the Y = 0 surface. 

Figure (3.1-1) now shows the grid in plane for a Cartesian coordinate 
system and illustrates the nomenclature described below. 



Fig. (3.1-1) Grid showing numbering, main control volumes and FORTRAN 
definition of grid quantities. (Cartesian coordinates) 
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The dashed lines in Figure (3.1-1) join the control -volume faces normal 
to the X and y directions. These are the control volumes for all 
dependent variables except the U and V velocity components. The control - 
volume faces pass mid-way between the grid nodes except near the boundaries 
where they pass through the boundary grid nodes. Thus, the control - 
volume faces always pass through points where the velocity component 
normal to the faces are stored. A normal velocity component at a control - 
volume face is presumed to prevail over that whole face. 

The control volumes for each of the velocity components U and V are 
displaced along the directions of these velocities. The control -volume 
faces normal to each of these directions pass through grid nodes on either 
side of the velocity component in question. Figure (3.1-2) illustrates 
this point. The geometric quantities associated with this grid system are 
defined in Table 3.1-1. 



V(I + I, J) 


Figure 3.1-2 


U- and V- Velocity Control Volumes 
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TABLE 3.1-1 


NO. 


GRID 

QUANTITY 


MEANING 


1 


X(I) 


Physical coordinate in the x-direction. 


2 


XDIF(I) 


The difference between X(I) and X(I-l); 
used as the distances s in calculating 
x-direction diffusion flux of <(>: r^A<()/6. 


3 


XS(I) 


The x-direction width of a main control 
volume. 


4 


XSU(I) 


The x-direction width of a U-velocity 
control volume. 


5 


Y(J) 


Physical coordinate in the y-direction. 


6 


YDIF(J) 


The difference between Y(J) and Y(J-l); 
used as the distances 6 in calculating 
y-direction diffusion flux of ((>. 


The y-direction width of a V-velocity 
control volume. 

The y-direction width of a V-velocity 
control volume. 




Because the faces of the main control volumes are stipulated to 
pass midway between the grid nodes, interpolation factors for the cal- 
culation of variables on these faces are equal to 0.5, except near 
boundaries where they are either 0 or 1. For the control volumes 
appropriate to U- and V- velocity components, however, the interpolation 
factors can differ from 0,5, if the spacing between grid nodes is chosen 
to be non-uniform. For this reason, interpolation factors are calculated 
in STRIDl and stored as FXP(I), FXM(I), FYP(J) and FYM(J). The sub- 
script refers to a grid node. The value of U for example at a grid node 
(I,J) is given by: 


FXP(I) * U(I + 1,J) + FXM(I) * U(I,J) (3.1-3) 


it is obvious from the above that FXM(I) is simply 1.0 - FXP(I), 

Calculation of the quantities tabulated above completes the tasks 
performed by STRIDl. STRIDl is called once from DAVE to calculate the 
initial grid quantities. Also, STRIDl is called from STRID3, after the 
W- velocity equation is solved and corrected for. 

3.2 Assembly of coefficients 

STRID2 

The STRID2 portion of STRIDA is used to calculate and store the 
convective mass velocities (i.e. pu, pv and pw in equation 1.2-1, 
crossing the control volume faces along the x and y directions). 

The arrays 6X and GY are used, respectively, to store these values. 

STRID2 is called, for each integration plane, once at the beginning 
of STRID3. Finally, STRID2 is called from STRID4, after the U- and V- 
velocities have been corrected. 

STRID3 

As referred to briefly above, STRID3 and STRIDA are the component 
subroutine of STRIDB which can be regarded as the main machinery of the 
HISS program. It is in STRID3, that the finite- difference equations 
appropriate to U,V,W and p' are assembled, in that order. After the 
appropriate equations are assembled, a call is made to SOMOD (a member 
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subroutine of ALLMOD) to determine whether any coefficients of the 
equations require modification. Once this is achieved, a call to 
SOLVE permits solutions of the finite-difference equations to be 
obtained. Once this preparation is completed, the coefficients of the 
pressure-correction equation set are assembled and the equations solved. 

STRID3 is then concluded with a section in which U's, Vs and p's are 
corrected. 

STRID4 

The assembly of and solutions to the finite-difference equations 
of all other dependent variables are performed in STRID4. The primary 
reason for keeping the functions of STRID3 and STRID4 apart is to permit 
iterations to be performed separately on each of these subprograms. The 
iterations themselves are initiated by the controlling subprogram MAIN. 

This iteration is used in HISS to obtain better distribution of K and e 
for the initial step. 

3. 3 Solution procedure for algebraic finite-difference equations 

SOLVE 

The function of the subprogram SOLVE is to arrange for the solutions 
to the finite-difference equations for each dependent variable NV, to be 
obtained. The solution procedure used is the application of a pair of Standard 
Tri-Diagonal Matrix (TPMA) traverses, one in each of the x and y directions. 

SOLVE has three major subdivisions. The first ends with the state- 
ment: 10 CONTINUE. It is in this part that the finite-difference coef- 
ficients are assembled in readiness for the subsequent operations. Fully- 
implicit practices are employed during the coefficient-assembly process. 

The second part of SOLVE is concerned with TDMA traverses in the x 
direction and ends with statement: 21 CONTINUE. 

The third and final part of SOLVE starts with this statement and 
concerns the TDMA traverses in the y direction. 

A call to SOLVE is made from STRID3 and STRID4, once for each dependent 
variable NV. This call achieves TDMA traverses in both the x and y direc- 
tions; however, which traverse is made first depends upon the value of 
the index IXY. The first traverse direction will be x if IXY = 1, and y 
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if IXY=2. IXY is set alternately to 1 or 2 by the statement: IXY=3-IXY 

which concludes both STRID3 and STRID4. The order of solution for problems 
for which HISS is applicable is not important. 

Along each traverse direction, the direction of sweep, i.e. 
whether from the first grid node to the last or vice versa, depends 
upon the value of an index, ISWP in the x direction and JSWP in the y 
direction. Each of these indices takes on a value 1 or 2 by statements 
which follow the end of parts 2 and 3 of SOLVE. For example, the state- 
ment following: 21 CONTINUE, reads: JSWP = 3 - JSWP. A value 1 implies 
that the sweep direction is from first to last grid node and 2 implies 
vice versa. 

On occasion, it may be required to perform more than one pair of 
TOMA traverses in a x-y plane, for any given dependent variable. The 
number of pairs of TOMA traverses is set by values assigned to an index 
array NSWP(NV). The program is set up with NSWP values equal to three 
except for the pressure-correction equation for which NSWP(NPP) = 6. 

3.4 Printout of field values of dependent variab les 

PRINTdSKIP, JSKIP) 

It is frequently required to print out the contents of the F array. 
This task is performed by subroutine PRINT. 

PRINT(ISKIP, JSKIP) provides for the printout of F(I,J,NV) where 
NV can attain a maximum value of NFPMAX, which is set as: 


NFPMAX = NNV + 3 


(3.4-1) 


The three extra values representing NVP, NRO and NMU, i.e. pressure, 
density and effective viscosity respectively. The decision as to whether 
a particular variable NV is printout or not, depends upon whether the 
corresponding IPRINT(NV) is equal to unity or not. 

The printout of each dependent variable NV is given a heading stored 
in TITLE( . . . ,NV). The formal parameters ISKIP and JSKIP permit the 
selective skipping of columns (I) and rows (J), when it is not required, 
for any reason, to printout the complete array of values of each variable. 
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4. PHYSICAL MODELLING 


4.0 Introduction 

This chapter describes the portion of HISS which calculates the 
physical properties required in the solution procedure. 

4.1 m 

The general policy is to confine all tasks associated with phy- 
sical modelling to subprogram AUX. Thus, the calculation of density p, 
effective diffusion coefficient r and sources and sinks s of the depen- 
dent variables are performed in separate member subroutines of AUX. 

There are five such subprograms in AUX, namely, DENSTY, VISCOS, GAMMA, 
SOURCE, and SPECIE. 

DENSTY 

This subroutine calculates densities pp u^RH0(I,J)j and p^^RH0D(I ,J)j 
in accordance with the discussion of section 1.7. 

VISCOS 

The function of VISCOS is to calculate the laminar and turbulent 
viscosity. In its present form, the major function is to store the 
turbulent viscosities in the array AMUT(J ,J ). . 

VISCOS also performs the function of calculating for the N and S 
boundaries, the effective boundary diffusion based upon the semi-logarith- 
mic law-of-the-wall . These diffusion coefficients are stored respectively 
in the arrays GAMN(I) and GAMS(I). 

GAMMA 

GAMMA is used to set values to the array GAM(I,J). GAM(IJ) is 
calculated from the formula appropriate to the turbulence model [4]. A 
call to GAMOD is then made in order to permit any modifications to be 
made to the GAM array. 

SOURCE 

Subroutine SOURCE is used to fill the arrays SU(I,J) and SP(I,J), 
there being a separate section for so doing for each dependent variable. 
The terms that fill these arrays are finite-difference equivalents to 
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the source/sink terms that appear in Tables 1.2-1. 

SPECIE 

Subroutine SPECIE is used to calculate the specie mass fraction 
from knowledge of the element mass fractions, pressure and temperature. 
Appendix A gives a full development of the technique employed. 
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5. PROBLEM DEPENDENT SECTIONS 


5.0 Introduction 

This chapter completes the detailed description of the functions 
of the various portions of HISS by discussing the subroutines which 
must be changed in order to specify a particular problem. 

5. 1 General policy 

As mentioned briefly earlier, the main machinery of the numerical 
calculation procedure is devoid of any problem-specification information. 
As a general policy, such information is provided through subroutines 
DATA and ALLMOD. 

5.2 BLOCK DATA 

BLOCK DATA serves to provide values to fluid properties, grid 
distributions, program control parameters and other information specific 
to each new problem, via DATA statements. The use of BLOCK DATA permits 
the program to be run with compilers common to both CDC and IBM machines. 

5.3 ALLMOD 

Subroutine ALLMOD is composed of five member subprograms BEGIN, 
GAMOD, SOMOD, UPSTRM and INJMOD. 

BEGIN 

The primary purpose of this subroutine is to provide initial values 
to all the dependent variable arrays, fluid-property arrays and other 
auxiliary arrays. The secondary purposes include provision of 'fixed' 
boundary conditions on the four boundaries of the calculation domain, and 
calculation of some auxiliary information required in the initializing 
process. 

GAMOD 

The function of modifying values of the array GAM(I,J) can be per- 
formed in GAMOD. Often, it is necessary to change only the boundary 
values of GAM; for instance, the diffusion coefficients relating to the 


37 



I I I 1 1 I 


II I I II III 


1 1 II 


velocity component normal to a wall boundary are set to zero at the 
boundary itself. The incorporation of wall functions, however, is not 
performed in GAMOD; this function being performed in SOMOD. 

SOMOD 

Any desired modifications to the finite-difference coefficients 
AXM(I,J), AXP{I,J) etc., and to the source terras SU(I,J) and SP(I,J), 
can be performed in SOMOD. There is a section in SOMOD corresponding 
to each dependent variable. As mentioned briefly earlier, the provision 
of wall-functions are made through SOMOD. This is achieved as follows: 
first, coefficients linking boundary grid nodes with their immediate 
neighbors inside the calculation domain are set to zero for each variable. 
Then if an index, for example KBCS for the South boundary, is set equal 
to unity, the appropriate wall flux is calculated and fed in through 
SU and SP. This calculation uses the value of GAMS(I), and the corresponding 
flux, FLUXS(I,NV) of the dependent variable is stored for purposes of 
printout. If the index is other than unity, no change is made to the 
coefficients. A similar indice, KBCN, is used for the ^orth wall. Such 
checks and modifications are made for the variables u,v,w,k,e,h, and f. 

UPSTRM 

The subprogram UPSTRM makes provisions to store upstream (i.e. 
previous integration plane) values of major variables. These include 
upstream members of the F array, stored as FU (I,J,NV), and certain other 
variables used in the calculations. Such a storage is required primarily 
for situations where iterations are resorted to and for calculation of 
source terms. 

INJMOD 

This subprogram is primarily used to specify the correct flux of a 
given variable through the array FINJ(NV). This array is used to modify 
the source term in SOMOD when injection occurs. 
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6. RESULTS AND CONCLUSIONS 


6.0 Introduction 

In this chapter comparison of predictions with experiment are given. 

further, sample results for two typical cases of interest in th.is contract 
are presented. 

6.1 Comparison to Experimental Data 

Figure 6.1-1 shows a comparison of results computed under this 
contract versus experimental data from NASA [7]. These results are 
for normal injection of hydrogen into a supersonic air stream. The 
conditions are defined in reference [7 ]. The main parameters are as 
follows : 


Air Stagnation Temperature = 

SOO'-K 

Air Stagnation Pressure 

1 . 38 MN/m2 

Air Mach Number = 

4.05 

Number of Jets = 

5 

Injector Diameter = 

.1 cm 

Jet Stagnation Pressure = 

.28 MN/m2 

Jet Stagnation Temperature = 

295°K 

Jet Mach Number = 

1 


The plot is made for a vertical plane aligned with the main flow direction 
and located at the center of the plate width. Three profiles of hydrogen 
concentration are shown at 30, 60 and 90 jet diameters downstream of the 
jet. Inspection of the figure shows excellent comparison between experi- 
mental and analytical results. The largest deviation occurs nearest the 
injector and probably can be ascribed to the elliptic nature of the flow 
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in the jet region. As discussed earlier the flow was forced to be para- 
bolic in this region. Thus, as the distance from the region increases, 
where the flow is parabolic, agreement between calculation and data is 
better. 

Also shown in Figure 6.1-1 is a plot of the experimental and computed 
shock location. The computed shock location was taken as occurring at the 
peak pressure at a given station between the wall and the outer edge of 
the computation region. Again this comparison shows excellent correlation 
with experiment. 

The correlation discussed above is a strong argument for the validity 
of the main results of this work which are described in the remaining 
part of this chapter. The cold flow case contains most of the physics of 
the cases of interest, namely supersonic compressible flow with a subsonic 
inner layer and strong cross stream interactions. 

6.2 Presentation of Typical Results 
6.2-1 Introduction 

Results for ten cases were calculated for this contract by the techni- 
que described in the first five chapters. For each case profiles of 19 
variables described earlier were computed for a grid of 12 grid lines 
transverse to the flow and 20 grid lines normal to the flow for a total 
of 240 points at each station in the main flow direction. Computations 
were made for several hundred stations. Thus, the total number of 
variables computed approaches nearly one-million for each case. In this 
section results for two typical cases will be presented. These results 
consist of hydrogen concentration, temperature and pressure distributions 
at three axial locations and at two locations transverse to the main flow 
direction. 
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6.2-2 Definition of Cases 

One case for normal injection (case 1) and one case for parallel 
injection (case 9) are described in this chapter. The geometrical con- 
figuration of these cases are shown in Figures 6.2-1 and 6.2-2 respectively. 
Table 6.2 defines the flow and thermal properties. 

Table 6.2 Definition of Properties 


Location 

Property (units) 

Case 1 

Case 9 

Main Stream 

Flow speed (m/s) 

675 

1585.2 

Condi tions 





Temperature (“K) 

70.6 

1178.6 


Mass fraction Ng 

.7676 

.487 


Mass fraction O 2 

.2325 

.263 


Mass fraction HgO 

0.0 

.25 


Pressure (N/m^) 

8720.0 

179300.0 

Jet 

Flow speed (m/s) 

1210.0 

2039.4 

Conditions 





Temperature (°K) 

250.0 

150.3 


i Mass fraction H 2 

1.0 

1.0 


Pressure N/m^ 

212000.0 

179300.0 


Mass flow kg/sec 

.000165 



.0109 


6.3 Graphs of Results 

All data are presented for planes normal to the wall and parallel with 
the main flow direction. One of these planes is located in the transverse 
direction directly at the hydrogen jet centerline and the other is located 
between jets (see Figures 6.2-1 and 6.2-2). The comparison of profiles at 
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Figure 6.2-2 Definition of Geometry and Planes of Profiles for Case 9. 




a given station in these two planes indicates the mixing rate in the trans- 
verse direction. Comparison of profiles in a given plane at various main 
flow stations indicates the rate of movement of hydrogen in a direction 
normal to the flow. 

Figures 6.3-1 through 6.3-3 are for normal injection (case 1). 

Figure 6.3-1 shows the distribution of hydrogen in any form at three axial 
locations {,118m, .216m and .246m). The first profile is before the injec- 
tion point (.196m) and therefore no hydrogen is present. At the .216m 
location, the hydrogen concentration levels (less than .01) indicate that 
substantial mixing has occurred, and the profiles indicate that mixing 
is primarily in the y direction. (The close spacing of the adjacent 
injectors limits x direction mixing.) Examination of the results at 
the .246m station shows continuing movement of hydrogen upward with an 
attendant smoothing of the profi'le, as would be expected. 

The temperature profiles of Figure 6.3-2 indicate that after injection 
and combustion that little transverse thermal gradient exists. This is 
partly due to the fact that the flame front is at the outer layer of the 
hydrogen zone which tends to diffuse the temperature in the transverse 
direction. Also it is seen that the reaction of hydrogen moves the thermal 
boundary layer outward from the plate with movement in the main flow 
direction. 

The pressure distribution for case 1 is shown in Figure 6.3-3. It 
exhibits the qualitative behavior one would expect. There is little 
transverse pressure difference due to the low velocities in that direction. 
A pressure spike can be seen downstream of the jet indicating a shock. 

The pressure below the shock is essentially uniform at a higher value 
than the free stream. It must be emphasized that the pressure distribu- 
tion in the jet region is subject to error due to the fact that the flow 
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Figure 6.3-3 Pressure Distribution for Case I 



(N/m*) 
I 


z = .246 m 



was forced to be parabolic when in actuality there is some upstream effect 
caused by the jet. 

Figures 6.3-4 through 6.3-6 are results for case 9 which is parallel 

injection. The diameter of the jets and jet spacing are much larger than 

for case 1. The profiles at injection (Z=0) show that there is little 

mixing of hydrogen between jets. However, as the flow moves forward the 
mixing increases such that the concentration of hydrogen between jets is 

roughly 50% of the hydrogen concentration in line with the jet at Z=.36 meter. 
However, these figures show that in the normal direction the movement of 
hydrogen upward between jets is not very effective. This is due to the low 
normal velocities for the parallel injection case. The temperature distri- 
bution in Figure 6.3-5 indicates a higher temperature between jets even 
though the hydrogen concentration is lower. This is because the combustion 
occurs at the fringe of the hydrogen zone. The combustion zone upper surface 
is indicated by the "spikes" in the temperature profile. This combustion is 
seen to grow with movement downstream. The pressure distribution shown in 
Figure 6.3-6 does not indicate any shocks. At the injection station (2=0), 
the flow expands in the low density region. However, at downstream locations 
the pressure becomes uniform at approximately the free stream value. 

6.4 Conclusions and Recommendations 

Qualitatively the results calculated by the method developed herein 
appear to be correct. Further comparison with cold flow data also show 
excellent agreement. The computational time for these results is excellent 
(1 minute on a CDC 7600). 

The assumption of parabolic flow in the injection region needs improve- 
ment because both data and calculations made in this work indicate that the 
pressure distribution is affected upstream by the jet. Recirculation, 
however, is probably not important. What is needed is a technique to allow 
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the pressure to be calculated via an elliptic procedure and the velocity by 
a parabolic one. Such a technique is advantageous over fully elliptic 
procedure because storage and computation time are much lower. 
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8. NOMENCLATURE 


E 

f 


h 


k 

K 



P 

^<i> 

P' 

Pr 

R 


T 


u, u* 

v, v* 

w, w* 


W 

X 

y 


turbulence constant (Table 1,2-2) 
specific heat 

turbulence constant (Table 1.2-2) 

turbulence constant (Table 1.2-2) 

law of wall constant (Table 1.2-2) 
mass fraction hydrogen in any form 

stagnation chemical enthalpy (Em^ + q 

2 p 

enthalpy of formation 
kinetic energy of turbulence 
constant in law of wall (Table 1.2-2) 
mass fraction of specie i 

pressure 

parameter defined by Equation (1.6-6) 

pressure correction 
Prandtl Number 
universal gas constant 
source term for 4> equation 
temperature 

corrected and uncorrected velocity in x direction 

corrected and uncorrected velocity in y direction 

corrected and uncorrected velocity in x direction 

molecular weight 

coordinate defined in Figure 1,2-1 
coordinate defined in Figure 1.2-1 
coordinate defined in Figure 1.2-1 
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NOMENCLATURE (Continued) 


conduction transport coefficient in (|) equation 

e dissipation rate of turbulence 

u laminar viscosity 

effective viscosity (laminar + turbulent) 

turbulent viscosity 

p density 

0, laminar Prandtl or Schmidt Number for * 

<P 

04. ^ turbulent Prandtl or Schmidt Number for d> 
t,(f> 

<fi general dependent variable in conservation equation 

SUBSCRIPT 

ref reference value 
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APPENDIX A 


A.O Background 

Reference [6] describes an equilibrium chemistry model developed and 
used in conjunction with a computer program to predict the properties in a 
hydrogen-oxygen flame. The main features of the model are described below. 
Four equilibrium reactions are assumed as follows: 


0 + 

0 

O 2 

(1) 

H + 

H 

H 2 

(2) 

0 + 

OH 

H 2 O 

(3) 

0 + 

H ^ 

OH 

(4) 


The six species involved in these reactions are considered to be present 
with nitrogen which is inert. In developing the equations to predict the 
equilibrium concentration of the species, two quantities are defined, 
namely 

Wn Wo 


F = mn^ + mn + mn^Q + mg^ (6) 

where X is the total fraction of oxygen in any form and F is the total 
fraction of hydrogen in any form. Since the molecular weight of the various 
oxygen species is approximately equal to that of nitrogen it is assumed that 
the rate of diffusion of nitrogen is equal to that of the oxygen; and, there- 
fore, nitrogen is present at any location in a fixed ratio to the fraction of 
oxygen compounds. This fraction, OFAC, is assumed constant and equal to the 
fraction of oxygen in tPie air being used as th,e oxidizer. Thus, 

X = (0FAC)(X + mM ) (7) 

‘V 

The total mass fraction must be unity which gives 


X + mw + F = 1 
2 


( 8 ) 
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Substituting equation (7) into equation (8) and solving for X gives 
X = (0FAC)(1 - F) (9) 

Therefore, if F is known, X can be determined by equation (9). 

From thermodynamic considerations the equilibrium constant is defined 
for the reaction aA + bB cC by 


Kp = (10) 

^ ya yb 
^A 

where P is in atmospheres. 

For each of the four reactions in the present model c-a-b = -1 
It is convenient to express the concentrations in terms of mass fractions. 
Noting that 

m^ = — Xj (11) 


and substituting equation (11) into equation (10) gives 



KpPW 


“I“b 



'"a"b 


( 12 ) 


Thus the equilibrium equations for the reactions (1 - 4) can be written 


K^ 



K2 



Ks 


"’H2O 

mo moH 


(13) 

(14) 

(15) 


= '^OH 

mo mj^ 


(16) 
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The condition of equilibrium is expressed by using four equilibrium 
constants for the four chemical reactions. Once F and X are specified, 
the six species (m[^^, m^, mg^, mo, mQH> q) can be determined from 
equations (5), (6), (13), (14), (15), and (16) in which these species 
are unknowns. The procedure adopted in [1] was to solve a differential 
equation for the distribution of F and then to solve equations (5), (6), 
(13), (14), (15), and (16) for the concentration of individual species. 

In reference [1] the values of species concentration just upstream of 
the point being considered were used to choose the largest term in both 
F and X. The mass fractions of remaining species were calculated from 
equations (13 - 16) using the upstream temperature for specifying the 
equilibrium constants and the values of the upstream species mass fraction 
which constitute the largest terms in F and X. New values of mass 
fraction of the largest species were calculated from equations (5) and (6) 
using the newly calculated values for the remaining species. The process 
was repeated until convergence to a specified limit was achieved. There 
are at least three major deficiencies with this approach: 

(1) Extreme care must be used in specifying the initial 
mass fractions; otherwise, the iteration procedure 
diverges; 

(2) Near the stoichiometric point several of the species have 
mass fractions which are approximately equal. This 
causes different species to have the largest mass fraction 
in successive steps in the iteration procedure resulting 
in instability; and 

(3) A large amount of computer logic is required to handle 
the extreme variations in the input parameters resulting 
in relatively large computational times. 

The continued use of the model and equations employed in reference [1] is 
recommended in this paper. However, a new solution procedure is proposed. 

The remaining discussion defines the proposed new solution procedure 
and gives some results and conclusions regarding its application. 
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A.l Deri vation of Solution Procedure 


The present procedure is based on the reduction of the number of 
variables under consideration. Two equations are derived as follows. 
Equations (13 - 16) are solved to give the relationships: 



(17) 


(18) 

Ka K4 

(19) 


(20) 

It is convenient to define the following parameters: 

A = 

(21) 

DOI 

(22) 

- _ K| K4 
Ki ^ K2 

(23) 

K3 

D = — 

Ki k; 

(24) 

Using equations (17 - 20) to eliminate mQ, 
equations (5) and (6) gives: 

rnn^O* '"H "’OH 
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( 25 ) 



where the definitions given by equations (21 - 24) have been used. 
These equations have the form of a quadratic equation: 


and solution [3] 


au^ + bu + c = 0 


u 


-2c 

B + - 4ac 


Note that a and 5 are always positive and c is always negative. Since 
u>0,the physically meaningful root is 


u 


-2c 

B + ^ b^ - 4ac 


(27) 


This particular form of quadratic expression is chosen since it does not 
require subtraction and gives greater precision. Using equation (27) to 
express the solution of equations (25) and (26) gives 




Equations (28) and (29) contain the unknowns and Before detailing 

the solution procedure several properties of these equations will be 
discussed. Table I gives the values of Kf, K 2 , K 3 , K 4 , A, B, C, and D for 
the temperatures from 100 to 6000 K at a constant PW product of 2.5 (cor- 
responding approximately to 1/10 atmosphere). Note that while the individual 
constants vary from less than 1 to 10 ^°°, the groups appearing in equations 
(28) and (29) vary much less and can be easily processed on a digital com- 
puter. Second it is to be observed that the equations have been derived so 

that a "natural" relationship has been established between mn and X, and 

2 

mH 2 and F, i.e. as X^, i" 02 '^» should also be observed 

that m 02 and nij^^ can never be less than zero or greater than 1. Finally, 
neither equation becomes indefinite as mn 2 or mo 2 approach zero. Thus, the 
equations are well-behaved and can be readily solved for wide variations in 
temperature and pressure. 

The solution procedure is described next. A value of mu is guessed 

2 

in the following way. If the value at the equivalent upstream station is 
known, then that value is used; otherwise, a value of zero will always 
lead to a converged solution. This assumed value of m^ is substituted 
into equation (28) which yields mn . Then the computed value of mn is 

'-'2 '“'0 

substituted into equation (29) which allows calculation of a new vaTue of 
m^^. The assumed and calculated values of mi^^ compared. If their 
values differ by more than a specified convergence criteria, the calculated 
value of m ^2 is taken as the assumed value and the process described above 
is repeated until convergence is obtained. The behavior of this solution 
technique is shown in Figure A. 1-1. 

Figure I shows with the broken line a plot of m^^ calculated versus 
m ^2 assumed. The correct value is achieved when the two values are equal. 

The locus of points for this situation is a straight line with a slope of 
unity. If m ^2 assumed is less than the correct value, the figure shows 
the calculated value will always be larger. Thus, when this calculated 
value of m ^2 ''s taken as the assumed value, the resulting newly calculated 
value will be closer to the correct one. The same argument can be made to 
show that if the initial assumed value of m^^ is too large, the iteration 
process will again cause convergence to the correct value. 
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A. 2 Typical Results 

The solution procedure described above was programmed on a digital 
computer and used to determine equilibrium concentrations for variation 
in F and X from zero to one for temperatures from 200 to 5700 K. The 
product PW was taken at 2.5 (approximately 1/10 atmosphere pressure). 

Table II shows these calculated results. The residuals in F and X are 
shown in the last two columns. Their values indicate that great precision 
can be obtained with this method. The convergence criteria for these 
calculations was 

l'"H2calculated " '^^^2®ssumedl ^ "’H2calculated 

Greater precision could easily be achieved by a stricter convergence 
criteria; however, the present one is sufficient for most practical cal- 
culations. 

It is concluded that the method offers a reliable, simple and extremely 
fast technique for solving the equilibrium equations arising from the 
chemical model considered in this work. 
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TABLE I 


Temp. Ki 


200 

5.0+100 

3. lE+99 

400 

2.8E+52 

2.7E+58 

600 

2.3E+33 

4. 3E+36 

800 

5.9E+23 

4.9E+25 

1000 

9.7E+17 

1.2E+19 

1200 

1.2E+14 

5.0E+14 

1400 

2.1E+11 

3.5E+11 

1600 

1.7E+09 

1.5E+09 

1800 

3.9E+07 

2.1E+07 

2000 

1.9E+06 

7.0E+05 

2200 

1.5E+05 

4.3E+04 

2400 

1.9E+04 

4.2E+03 

2600 

3.4E+03 

5.8E+02 

2800 

7.5E+02 

l.OE+02 

3000 

2.0E+02 

2.4E+01 

3200 

6.3E+01 

6.8E+00 

3400 

2.3E+01 

2.1E+00 

3600 

9.3E+00 

7.9E-01 

3800 

4.1E+00 

3.2E-01 

4000 

1.9E+00 

1.4E-01 

4200 

l.OE+00 

6.8E-02 

4400 

5.5E-01 

3.5E-02 


K3 

Ku 

A 

2.6+100 

2.6+100 

4.4E-51 

3.2E+51 

6.2E+59 

5.9E-27 

5.8E+32 

9.0E+37 

2.0E-17 

2.2E+23 

9.3E+26 

1.2E-12 

4.5E+17 

2.2E+20 

l.OE-09 

7.0E+13 

8.2E+15 

8.7E-08 

1.2E+11 

5.4E+12 

2.1E-06 

l.lE+09 

2.2E+10 

2.4E-05 

2.7E+07 

3.0E+08 

1.5E-04 

1.4E+06 

9.8E+06 

7.2E-04 

1.2E+05 

5.8E+05 

2.5E-03 

1.6E+04 

5.5E+04 

7.0E-03 

2.8E+03 

7.6E+03 

1.7E-02 

6.4E+02 

1.3E+03 

3.6E-02 

1.7E+02 

3.1E+02 

7.0E-02 

5.7E+01 

8. 5E+01 

1.2E-01 

2.1E+01 

2.7E+01 

2.0E-01 

8.6E+00 

9.8E+00 

3.2E-01 

3.8E+00 

3.9E+00 

4.9E-01 

1.8E+00 

1.7E+00 

7.0E-01 

9.8E-01 

8.3E-01 

9.8E-01 

5.4E-01 

4.2E-01 

1.3E+00 


B' 

C 

D 

1.7E-50 

2.1E+00 

2.5E+50 

6.0E-30 

l.lE-04 

4.3E+29 

4.7E-19 

5.7E-03 

l.OE+19 

1.4E-13 

4.0E-02 

4.8E+13 

2.7E-10 

1.2E-01 

2.8E+10 

4.4E-08 

2.7E-01 

1.9E+08 

1.6E-06 

4.7E-01 

5.5E+06 

2.5E-05 

7.0E-01 

3.8E+05 

2.1E-04 

9.5E-01 

4.6E+04 

l.lE-03 

1.2E+00 

8.7E+03 

4.8E-03 

1.4E+00 

2.2E+03 

1.5E-02 

1.7E+00 

6.9E+02 

4.1E-02 

2.0E+00 

2.6E+02 

9.6E-02 

2.2E+00 

l.lE+02 

2.0E-01 

2.5E+00 

5.5E+01 

3.8E-01 

2.7E+00 

2.9E+01 

6.7E-01 

2.9E+00 

1.6E+01 

l.lE+00 

3.1E+00 

l.OE+01 

1.7E+00 

3.3E+00 

6.5E+00 

2.6E+00 

3.5E+00 

4.4E+00 

3.8E+00 

3.7E+00 

3.0E+00 

5.3E+00 

3.8E+00 

2.2E+00 


65 


TABLE I. 
(Continued) 


Temp. 

Ki 

K2 

4600 

3.2E-01 

1.8E-02 

4800 

1.9E-01 

l.OE-02 

5000 

1.2E-01 

6.4E-03 

5200 

7.9E-02 

4.0E-03 

5400 

5.3E-02 

2.5E-03 

5600 

3.6E-02 

1.7E-03 

5800 

2.6E-02 

l.lE-03 

6000 

1.8E-02 

8.1E-04 


Ka K4 a 


3.1E-01 

2.2E-01 

1.7E+00 

1.9E-01 

1.3E-01 

2.2E+00 

1.2E-01 

7.7E-02 

2.8E+00 

7.8E-02 

4.8E-02 

3.5E+00 

5.3E-02 

3.1E-02 

4.3E+00 

3.6E-02 

2.0E-02 

5.2E+00 

2.6E-02 

1.4E-02 

6.1E+00 

1.8E-02 

9.8E-03 

7.2E+00 


B C D 


7.2E+00 

4.0E+00 

1.6E+00 

9.6E+00 

4.1E+00 

1.2E+00 

1.2E+01 

4.2E+00 

9.5E-01 

1.5E+01 

4.4E+00 

7.5E-01 

1.9E+01 

4.5E+00 

6.0E-01 

2.4E+01 

4.6E+00 

4.9E-01 

2.9E+01 

4.7E+00 

4.1E-01 

3.5E+01 

4.8E+00 

3.4E-01 
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TABLE II 


Temp . 

F 

X 

CM 


■"h 

mo 

■"oh 

"’ H20 

RES * ** 

RE$2 

200 

0 

.23200 

.00000 

.23200 

0 

0 

0 

.00000 

-.00000 

0 

200 

.2 

. 18560 

. 17677 

.00000 

0 

0 

0 

.20903 

-.00020 

0 

200 

.4 

.13920 

.38259 

.00000 

0 

0 

0 

. 15665 

-.00005 

0 

200 

.6 

.09280 

.58840 

.00000 

0 

0 

0 

. 10441 

-.00001 

0 

200 

.8 

.04640 

.79419 

.00000 

0 

0 

0 

.05228 

-.00008 

0 

200 

1.0 

.00000 

1.00000 

.00000 

0 

0 

0 

.00001 

-.00000 

0 

700 

0 

.23200 

.00000 

.23200 

0 

0 

0 

.00000 

-.00000 

0 

700 

.2 

. 18560 

.17677 

.00000 

0 

0 

0 

.20903 

-.00020 

0 

700 

.4 

.13920 

.38259 

.00000 

0 

0 

0 

.15665 

-.00005 

0 

700 

.6 

.09280 

.58840 

.00000 

0 

0 

0 

.10441 

-.00001 

0 

700 

.8 

.04640 

.79419 

.00000 

0 

0 

0 

.05228 

-.00008 

0 

700 

1.0 

.00000 

1.00000 

.00000 

0 

0 

0 

.00001 

-.00000 

0 

1200 

0 

. 23200 

.00000 

.23200 

0 

0 

0 

.00000 

-.00000 

0 

1200 

.2 

.18560 

.17677 

.00000 

0 

0 

0 

.20903 

-.00020 

0 

1200 

.4 

.13920 

. 38259 

.00000 

0 

0 

0 

. 15665 

-.00005 

0 

1200 

.6 

.09280 

.58840 

.00000 

0 

0 

0 

.10441 

-.00001 

0 

1200 

.8 

.04640 

.79419 

.00000 

0 

0 

0 

.05228 

-.00008 

0 

1200 

1.0 

.00000 

1.00000 

.00000 

0 

0 

0 

.00001 

-.00000 

0 

1700 

0 

.23200 

.00000 

.23196 

.00000 

.00004 

0 

.00000 

-.00000 

0 

1700 

.2 

.18560 

.17675 

.00000 

.00003 

.00000 

0 

.20902 

-.00020 

0 

1700 

.4 

.13920 

.38255 

.00000 

.00004 

.00000 

0 

.15665 

-.00005 

0 


* Residual in F equation (29 

** Residual in X equation (28 
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TABLE II . 
( Continued ) 


Temp . 

F 

X 


% 

niH 

nio 

"oh 

" H2O 

RES * 

1 

RES ** 

2 

1700 

.6 

.09280 

.58835 

.00000 

.00005 

.00000 

0 

.10441 

-.00001 

0 

1700 

.8 

.04640 

.79413 

.00000 

.00006 

.00000 

0 

.05228 

-.00008 

0 

1700 

1.0 

.00000 

.99993 

.00000 

.00007 

.00000 

0 

.00000 

- . 00000 

0 

2200 

0 

.23200 

.00000 

.22970 

.00000 

,00230 

.00000 

.00000 

-.00000 

0 

2200 

.2 

. 18560 

.17574 

.00000 

.00105 

.00000 

.00034 

.20867 

-.00020 

0 

2200 

.4 

. 13920 

.38106 

.00000 

.00155 

.00000 

.00017 

.15647 

-.00005 

0 

2200 

.6 

.09280 

. 58648 

.00000 

.00192 

.00000 

.00009 

.10431 

-.00001 

0 

2200 

.8 

.04640 

.79196 

.00000 

.00223 

.00000 

.00004 

.05224 

-.00008 

0 

2200 

1,0 

.00000 

.99749 

.00000 

.00251 

.00000 

.00000 

.00000 

-.00000 

0 

2700 

0 

. 23200 

.00000 

.20310 

.00000 

.02890 

.00000 

.00000 

-.00000 

0 

2700 

.2 

. 18560 

. 16687 

.00005 

.01034 

.00045 

.00620 

.20189 

-.00019 

0 

2700 

.4 

. 13920 

.36746 

.00001 

.01534 

.00016 

.00317 

.15311 

-.00004 

0 

2700 

.6 

.09280 

. 56941 

.00000 

.01910 

.00007 

.00170 

.10253 

-.00001 

0 

2700 

.8 

. 04640 

.77200 

.00000 

.02223 

.00002 

.00073 

.05148 

-.00007 

0 

2700 

1.0 

.00000 

.97501 

. 00000 

.02499 

.00000 

.00000 

.00000 

-.00000 

0 

3200 

0 

.23200 

.00000 

.10689 

.00000 

.12511 

.00000 

,00000 

-.00000 

0 

3200 

.2 

. 18560 

.13505 

.00142 

.04597 

. 01444 

.03815 

.15067 

-.00011 

0 

3200 

.4 

.13920 

.31439 

.00019 

.07014 

.00527 

.02124 

.12800 

-.00003 

0 

3200 

.6 

.09280 

. 50086 

.00004 

.08853 

.00231 

.01175 

.08933 

-.00001 

0 

3200 

.8 

.04640 

.69065 

.00001 

.10395 

.00086 

.00514 

.04586 

-.00006 

0 

3200 

1.0 

.00000 

. 88249 

.00000 

.11751 

.00000 

.00000 

.00000 

-.00000 

0 


* Residual in F equation ( 29 ) 

** Residual in X equation ( 28 ) 
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TABLE II. 
(Continued) 


Temp. 

F 

X 



mn 

mo 

'"oh 

'"H2O 

RES* 

RES** 

2 

3700 

0 

. 23200 

.00000 

.02210 

.00000 

.20990 

.00000 

.00000 

-.00000 

0 

3700 

.2 

. 18560 

.07874 

.00406 

.11326 

.08998 

.05854 

.04102 

-.00000 

0 

3700 

.4 

.13920 

.20737 

.00102 

.18380 

.04514 

.04766 

.05420 

-.00000 

0 

3700 

.6 

.09280 

.35319 

.00025 

. 23987 

.02252 

.03103 

.04604 

-.00010 

0 

3700 

.8 

.04640 

.50841 

.00004 

.28779 

.00897 ■ 

.01483 

.02640 

-.00002 

0 

3700 

1.0 

.00000 

.66970 

.00000 

.33030 

.00000 

.00000 

.00000 

-.00000 

0 

4200 

0 

. 23200 

.00000 

.00357 

.00000 

.22843 

.00000 

. 00000 

-.00000 

0 

4200 

.2 

. 18560 

.02919 

.00167 

. 16888 

.15622 

.02599 

.00366 

-.00000 

0 

4200 

.4 

.13920 

.09410 

.00072 

.30323 

.10272 

.03068 

.00775 

-.00001 

0 

4200 

.6 

.09280 

.35319 

.00025 

.23987 

.02252 

.03103 

.04604 

-.00010 

0 

4200 

.8 

.04640 

.27764 

.00005 

. 52084 

.02758 

.01415 

.00614 

-.00000 

0 

4200 

1.0 

.00000 

.38593 

.00000 

.61407 

.00000 

.00000 

.00000 

-.00000 

0 

4700 

0 

.23200 

.00000 

.00076 

.00000 

.23124 

.00000 

.00000 

0 

0 

4700 

.2 

. 18560 

.00904 

.00045 

.19045 

.17716 

.00824 

.00027 

0 

0 

4700 

.4 

.13920 

.03336 

.00023 

. 36589 

.12760 

.01140 

.00072 

0 

0 

4700 

.6 

.09280 

.06985 

.00010 

.52942 

.08188 

.01059 

.00097 

0 

0 

4700 

.8 

.04640 

.11632 

.00002 

.68320 

.03949 

.00659 

.00078 

0 

0 

4700 

1.0 

.00000 

.17119 

.00000 

.82881 

.00000 

.00000 

.00000 

0 

0 

5200 

0 

. 23200 

.00000 

.00022 

.00000 

.23178 

.00000 

.00000 

0 

0 

5200 

.2 

.18560 

.00308 

.00013 

.19675 

. 18278 

.00283 

.00003 

0 

0 

5200 

.4 

.13920 

.01195 

.00007 

. 38780 

.13517 

.00413 

.00008 

0 

0 


* Residual in F equation (29) 

** Residual in X equation (28) 
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Temp 


5200 

5200 

5200 

5700 

5700 

5700 

5700 

5700 

5700 


* 

** 


TABLE II. 

(Continued) 

F X m^^ mQ^ m^ mg mpj^ RESi RE$2 


.6 

.09280 

.02615 

.00003 

. 57360 

.8 

.04640 

.04526 

.00001 

.75458 

1.0 

.00000 

.06891 

.00000 

.93109 

0 

.23200 

.00000 

.00008 

.00000 

.2 

. 18560 

.00122 

.00005 

. 19871 

.4 

. 13920 

.00483 

.00003 

.39507 

.6 

.09280 

.01075 

.00001 

.58916 

.8 

.04640 

.01889 

.00000 

.78105 

1.0 

.00000 

.02918 

.00000 

.97082 


.08889 

.00402 

.00011 

0 

0 

.04385 

.00261 

.00009 

0 

0 

.00000 

.00000 

.00000 

0 

0 

.23192 

.00000 

.00000 

0 

0 

. 18448 

.00113 

.00000 

0 

0 

. 13758 

.00168 

.00001 

0 

0 

.09121 

.00166 

.00002 

0 

0 

.04535 

.00109 

.00001 

0 

0 

.00000 

.00000 

.00000 

0 

0 


Residual in F equation (29) 
Residual in X equation (28) 
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APPENDIX B 


DEFINITION OF TERMS APPEARING IN COMMON STATEMENTS 


AGEOM 

AK 

AKFAC 

AMU 

AMUREF 

AMUT 

BXE 

BYN 

CD 

CDQR 

CDRT 

CDTQ 

CP 

CPDCV 

Cl 

C2 

DEN 

DZ 

EE 

EX 

FLOINJ 

FLUXN, FLUXS 

FRA 

FRAM 

FXM, FXP, FYM, FYP 
GAMN, GAMS 
GASCON 


factor in grid expansion defined in section 3.1 
mixing length constant, K 

factor relating turbulence energy to mean motion 
energy 

laminar viscosity, u 

reference laminar viscosity 

turbulent viscosity, 

width of flow region along x-coordinate 

width of flow region along y-coordinate 

constant in turbulence model, Cp 

Cq raised to 1/4 power 

Cq raised to 1/2 power 

Cq raised to 3/4 power 

specific heat 

specific heat at constant pressure divided by 
specific heat at constant volume 

constant in turbulence model , Cj 

constant in turbulence model, C 2 

reference density 

step length in z direction 

constant in law-of-the-wal 1 , E 

factor by which forward step size is increased 

mass flow rate injected by a jet 

flux of dependent variable at N and S boundaries 
respecti vely 

fraction of boundary height used in calculating 
forward step 

maximum fraction of boundary layer height used in 
calculating forward step 

interpolation factors to indicate distance of a 
pressure node from the neighboring velocity node 

boundary value of GAM at N and S boundaries, 
respecti vely 

universal gas constant 
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GREAT 

H 

HO 
IINJ 
I MAX 

INJSTP 

INTEX 

IPRINT 

I REF 
ISOLVE 

ISTEP 

ISTR 

ISWP 

ITERF 

ITERM 

IXY 

JM 

JMAX 

OSTR 

JSWP 

KBCN, KBCS 
L 

LASTEP 

LCV 

LPl 

M 

MCV 

MPl 

NFM 


large number used in limiting overflows, etc. 

recovery factor 

enthalpy of formation 

index to control entry to INJMOD 

maximum value of I for which storage locations 
are provided 

value of ISTEP at injection location 

index: 1 for internal flows: 2 for external flows 

values of a variable printed only if corresponding 
value of IPRINT is unity 

reference value of I for printout 

equation for a variable is solved only if corre- 
sponding value of ISOLVE is unity 

number of station 

value of I for the first storage location for a 
given variable 

index controlling direction if TDMA sweep along I 
coordi nate 

counter for iterations on the <() equation 

counter for iteration on the momentum equation 

index controlling the first direction of a TDMA 
sweep 

(J - DJMAX 

maximum value of J for which internal storage is 
provided 

value of J denoting the first internal storage 
location for a given variable 

index denoting the direction of sweep while per- 
forming TDMA sweep in y-di recti on 

index denoting nature of N and S boundaries, respec- 
tively: 1 = wall; 2 = symmetry; 3 = free. 

number of grid lines minus 1 in x-di recti on 

number of station for which integration is to end 

number of main control volumes in x-di recti on 

number of grid lines in x-di recti on 

number of grid lines minus 1 in y-di recti on 

number of main control volumes in y-di recti on 

number of grid lines in y-di recti on 

number of variable (NV)*IMAX*JMAX 


72 



NFPMAX 

NH 

NH2 

NH20 

MMU 

NNV 

NN2 

NO 

NOH 

N02 

NPJUMP 

NPP 

NRO 

NSWP 

NV 

NVD 

NVF 

NVH 

NVK 

NVP 

NVT 

NVU 

NVV 

NVW 

PIN 

PJAY 

PR 

PRLAM 

RETRAN 

SMALL 

TIN 

TINJ 


maximum number of variables for which printout 
is arranged 

identifier of atomic hydrogen concentration 

identifier of molecular hydrogen concentration 

identifier of water vapor concentration 

identifier of viscosity 

number of dependent variables 

identifier of molecular nitrogen concentration 

identifier of atomic oxygen concentration 

identifier of OH concentration 

identifier of molecular oxygen concentration 

value equals ISTEP for printout of profile of 
variables 

identifier of pressure correction 

identifier of density 

number of TOMA sweeps on a variable 

serial number of any variable 

identifier for dissipation rate of turbulence 

identifier for hydrogen concentration in any form 

identifier for enthalpy 

identifier for turbulent kinetic energy 

identifier for pressure 

identifier for temperature 

identifier for velocity in x-di recti on 

identifier for velocity in y-di recti on 

identifier for velocity in z-di recti on 

inlet pressure 

resistance of laminar sublayer, P, 

<p 

Prandtl number, Pr 

laminar Prandtl number/Schmidt number, 

transition Reynolds number 

small number used to control underflows, etc. 

inlet free stream temperature 

inlet jet temperature 




TITLE 

TWAL 

UIN 

VIN 

VINJ 

WIN 

WM 

X 

XDIF 

XS 

XSU 

Y 

YDIF 

YS 

YSV 

ZD 

ZETA 

ZINJ 

ZLAST 

ZRE 

ZU 


any title to describe calculations 

wall temperature 

inlet x-di recti on velocity 

inlet y-di recti on velocity 

velocity of jet injection 

inlet free-stream velocity in z-di recti on 

molecular weight 

x-coordi nates 

distance between grid lines in x-di recti on 
width of main control volume in x-di recti on 
width of u velocity control volume in x-di recti on 
y- coordinates 

distance between grid lines in y-di recti on 

width of main control volume in y-direction 

width of V velocity control volume in y-direction 

downstream location on control volume face in z- 
di recti on 

x/BXE 

z location of injection 

z location at end of integration 

values of z for which printout is desired 

upstream location of control volume face in z- 
di recti on 
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